Componentwise accurate Brownian motion computations using Cyclic Reduction

Markov-modulated Brownian motion is a popular tool to model continuous-time phenomena in a stochastic context. The main quantity of interest is the invariant density, which satisfies a differential equation associated with the quadratic matrix polynomial P (z) = V z−Dz+Q, where the matrices V and D are diagonal and Q is the transition matrix of a discrete-time Markov chain. Its solution is typically constructed by computing an invariant pair of P (z) associated with its eigenvalues in the left half-plane, or by solving the matrix equation XV − XD + Q = 0. We show that these tasks can be solved using a componentwise accurate algorithm based on Cyclic Reduction, generalizing the recently appeared algorithms for the linear case (V = 0). We give a proof of the numerical stability of our algorithm in the componentwise sense; the same proof applies to Cyclic Reduction in a more general M-matrix setting which appears in other applications such as the modelling of QBD processes.